rm(list = ls())
library(tidyverse)
library(lubridate)
library(ggplot2)
library(ggthemes)
library(estimatr)

df_cjt = readRDS("data/baseline_cjt.rds")

# Function ----------------------------------------------------------------
effect_plotter <- function(estimate_vec,
                           se_vec,
                           names.variables,
                           names.levels,
                           effect.label,
                           x.lower = NULL,x.upper = NULL, labs = F, title = "ABC"){
  
  if(length(names.variables) != length(names.levels)){
    stop("Number of sets of levels does not match number of variables!")
  }
  
  n.vars <- length(names.variables)
  name <- c()
  code <- c()
  group <- c()
  pe <- c()
  se <- c()
  solo <- c()
  
  for (i in 1:n.vars){
    
    nv <- names.variables[i]
    nl <- names.levels[[i]]
    
    if(length(nl) > 1){
      name <- c(name,nv,nl,NA)
    } else {
      name <- c(name,nv,NA)
    }
    
    code.tmp <- c()
    solo.tmp <- c()
    if (length(nl) > 1){
      code.tmp <- c(0,0,rep(1,length(nl) - 1),0)
      solo.tmp <- c(0,0,rep(0,length(nl) - 1),0)
    } else {
      code.tmp <- c(1,0)
      solo.tmp <- c(1,0)
    }
    code <- c(code,code.tmp)
    solo <- c(solo,solo.tmp)
    
    group.tmp <- c()
    if (length(nl) > 1){
      group.tmp <- c("empty",rep(paste("group",i,sep=""),
                                 length(nl)),"empty")
    } else {
      group.tmp <- c(1,"empty")
    }
    group <- c(group,group.tmp)
    
    pe.tmp <- c()
    se.tmp <- c()
    if (length(nl) > 1){
      pe.tmp <- c(NA,rep(0,length(nl)),NA)
      se.tmp <- c(NA,rep(0,length(nl)),NA)
    } else {
      pe.tmp <- c(0,NA)
      se.tmp <- c(0,NA)
    }
    pe <- c(pe,pe.tmp)
    se <- c(se,se.tmp)
    
  }
  
  code[code==1] <- seq(1:sum(code))
  pdat <- data.frame(name,group,code,pe,se,solo)
  pdat <- pdat[-nrow(pdat),]
  pdat$order <- rev(seq(1:nrow(pdat)))
  
  #Extract pe and se estimates
  pevec <- estimate_vec
  sevec <- se_vec
  
  #Input the estimates into the framework dataframe
  for (i in 1:max(pdat$code)){
    pdat$pe[pdat$code == i] <- pevec[i+1]
    # +1 because of intercept
    pdat$se[pdat$code == i] <- sevec[i+1]
  }
  
  pdat$name <- as.character(pdat$name)
  pdat$group <- as.character(pdat$group)
  pdat$name[is.na(pdat$name)] <- ""
  
  if (!is.null(x.lower) & !is.null(x.upper)){
    theplot <- plotit(d = pdat, effect.label = effect.label,
                      x.lower = x.lower, x.upper = x.upper, labs = labs, title = title)
  } else {
    theplot <- plotit(d = pdat, effect.label = effect.label, labs = labs, title = title)
  }
  
  return(theplot)
  
}



plotit <- function(d,effect.label,x.lower = NULL,x.upper = NULL, labs, title, facet = NULL){  
  
  if(labs == F){
    lab = "element_blank()"
  }else{
    lab = "element_text(size = base_size, hjust = 0 , vjust=.5, face = rev(plot.face))"
  }
  
  
  CIs <- function(d){
    d$upper <-d$pe + 1.96*d$se
    d$lower <-d$pe - 1.96*d$se
    return(d)
  }
  d<- CIs(d)
  
  plot.labels <- as.character(d$name)
  plot.labels[d$group != "empty" & 
                !is.na(plot.labels) & d$solo == 0] <- 
    paste("    ",
          plot.labels[d$group != "empty" & 
                        !is.na(plot.labels) & d$solo == 0], sep = "")
  plot.face <- ifelse(is.na(d$pe) | d$solo == 1, "bold", "plain")
  d$ref <- ifelse(d$pe != 0 | is.na(d$pe), "normal", "reference")
  
  d$order <- as.factor(d$order)
  
  #Construct plot
  p <- ggplot(d,aes(y=pe,x=order,color=group,size=0.3,shape=ref)) + 
    scale_shape_manual(values = c(16,1))
  
  if (!is.null(x.lower) & !is.null(x.upper)){
    p <- p + coord_flip(ylim = c(x.lower,x.upper)) 
  } else {
    p <- p + coord_flip() 
  }
  
  if(!is.null(facet)){
    p <- p + facet_wrap(facet)
  }
  
  
  p <- p + ylab(effect.label)
  p <- p + geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype=1) 
  p <- p + geom_pointrange(aes(ymin=lower,ymax=upper,width=.4),
                           position="dodge",size=.5)
  p <- p + scale_colour_discrete("Attribute:") + 
    scale_x_discrete(name="",labels=rev(plot.labels)) 
  
  theme_bw1 <- function(base_size = 6, base_family = "") {
    theme_economist(base_size = base_size, base_family = base_family)
    theme(
      axis.text.x = element_text(size = base_size, hjust = .5 , vjust=1, colour = "black"),
      axis.text.y = eval(parse(text = lab)),
      axis.ticks.y = element_line(colour = "grey50"),
      axis.ticks.x = element_blank(),
      axis.title.y = element_text(size = base_size+2,angle=90,
                                  vjust=.01,hjust=.1),
      axis.title.x = element_text(size = base_size),
      legend.position = "none"
    )
  }
  base_size = 8
  base_family = ""
  p <- p + theme_economist() + theme(axis.text.x = element_text(size = base_size, hjust = .5 , vjust=1, colour = "black"),
                                     axis.text.y = eval(parse(text = lab)),
                                     #axis.ticks.y = element_line(colour = "grey50"),
                                     axis.ticks.x = element_blank(),
                                     axis.title.y = element_text(size = base_size,angle=90,
                                                                 vjust=.01,hjust=.1),
                                     axis.title.x = element_text(size = base_size),
                                     legend.position = "none",
                                     plot.title = element_text(size = base_size+2, face = "bold", hjust = .5, vjust = 2)) #hjust centers title
  p <- p + theme(panel.grid.major = element_line(size = 0.28)) + 
    theme(panel.grid.minor = element_blank())
  p <- p + labs(title = title)
  theme(title = element_text(size = 5))
  
  return(p)
  
}

# Preparing for analysis --------------------------------------------------
df_cjt$cjt_services_lebanon %>% table

df_cjt = df_cjt %>% 
  mutate(cjt_safety2 = recode(cjt_safety, "All of Syria is safe" = "Hometown or Syria safe", "Your hometown is safe" = "Hometown or Syria safe"),
         cjt_safety2 = factor(cjt_safety2, levels = c("Your hometown is not safe", "Hometown or Syria safe")),
         cjt_safety = factor(cjt_safety, levels = c("Your hometown is not safe", "Your hometown is safe", "All of Syria is safe")),
         cjt_econ_syria_jobs = factor(cjt_econ_syria, levels = c("There are few job opportunities", "There are many job opportunities", "Public services hard to attain", "Public services easy to attain")),
         cjt_econ_syria_services = factor(cjt_econ_syria, levels = c("Public services hard to attain", "Public services easy to attain", "There are few job opportunities", "There are many job opportunities")),
         cjt_conscription = factor(cjt_conscription, levels = c("Military conscription remains", "Military conscription ended")),
         cjt_econ_lebanon_jobs = factor(cjt_econ_lebanon, levels = c("No good job in Lebanon", "Good job in Lebanon", "Public services not available in Lebanon", "Public services available in Lebanon")),
         cjt_econ_lebanon_services = factor(cjt_econ_lebanon, levels = c("Public services not available in Lebanon", "Public services available in Lebanon", "No good job in Lebanon", "Good job in Lebanon")),
         cjt_networks = factor(cjt_networks, levels = c("Friends elsewhere", "Friends in Lebanon", "Friends in Syria"))
  )

cjt_safety2 = unique(df_cjt$cjt_safety2)
cjt_safety = unique(df_cjt$cjt_safety)
cjt_econ_syria_jobs = unique(df_cjt$cjt_econ_syria)
cjt_econ_syria_services = unique(df_cjt$cjt_econ_syria)
cjt_conscription = unique(df_cjt$cjt_conscription)
cjt_econ_lebanon_jobs = unique(df_cjt$cjt_econ_lebanon)
cjt_econ_lebanon_services = unique(df_cjt$cjt_econ_lebanon)
cjt_networks = unique(df_cjt$cjt_networks)

refcats = tibble(variable = c("Your hometown is not safe", "There are few job opportunities", "Public services hard to attain", 
                              "Military conscription remains", "No good job in Lebanon", "Public services not available in Lebanon",
                              "Friends elsewhere"),
                 coef = 0, se = 0, attribute = 0)

attributes = tibble(variable = c("Safety", "Jobs in Syria", "Public Services in Syria", "Conscription", "Jobs in Lebanon",
                                 "Public services in Lebanon", "Networks"),
                    coef = NA_real_,
                    se = NA_real_, attribute = 1)

conjoint = df_cjt %>% filter(complete.cases(outcome, cjt_safety, cjt_econ_syria_jobs, cjt_econ_syria_services, cjt_conscription, cjt_econ_lebanon_jobs, cjt_econ_lebanon_services, cjt_networks))

basemod1 = lm_robust(outcome ~  cjt_safety + cjt_econ_syria_jobs + cjt_conscription + cjt_econ_lebanon_jobs + cjt_networks, weights = weights, clusters = response_num, data = conjoint) %>% tidy
basemod2 = lm_robust(outcome ~  cjt_safety + cjt_econ_syria_services + cjt_conscription + cjt_econ_lebanon_services + cjt_networks, weights = weights, clusters = response_num, data = conjoint) %>% tidy

pevec = c(basemod1[1:4,2], basemod2[4,2], basemod1[7:8,2], basemod2[8,2], basemod1[11:12,2])
sevec = c(basemod1[1:4,3], basemod2[4,3], basemod1[7:8,3], basemod2[8,3], basemod1[11:12,3])
pvalvec = c(basemod1[1:4,5], basemod2[4,5], basemod1[7:8,5], basemod2[8,5], basemod1[11:12,5])

main_reg = lm_robust(outcome ~  cjt_safety + cjt_econ_syria_jobs + cjt_conscription + cjt_econ_lebanon_jobs + cjt_networks, weights = weights, clusters = response_num, data = conjoint)

library(texreg)

effect_plotter(estimate_vec = pevec, se_vec = sevec, names.variables = c("Safety", "Jobs in Syria", "Public Services in Syria", "Conscription", "Jobs in Lebanon",
                                                                         "Public services in Lebanon", "Networks"), 
               names.levels = list(c("Hometown not safe", "Hometown safe", "All Syria safe"), 
                                   c("Few job opportunities", "Many job opportunities"),
                                   c("Unavailable public services", "Available public services"),
                                   c("Military conscription remains", "Military conscription ended"),
                                   c("Lack good job", "Possess good job"),
                                   c("Public services not available", "Public services available"),
                                   c("Friends/family elsewhere", "Friends/family in Lebanon", "Friends/family in Syria")
               ),
               effect.label = "Effect on probability of intending to return", 
               x.lower = -0.1, x.upper = 0.45, labs = T, title = "")

ggsave("figures/conjoint.png", height = 5, width = 5)




# By Safety ---------------------------------------------------------------

conjoint = df_cjt %>% filter(complete.cases(outcome, cjt_safety2, cjt_econ_syria_jobs, cjt_econ_syria_services, cjt_conscription, cjt_econ_lebanon_jobs, cjt_econ_lebanon_services, cjt_networks))

conjoint %>% filter(cjt_safety2 == "Your hometown is not safe" & cjt_conscription == "Military conscription remains") %>% pull(outcome) %>% table %>% prop.table

conjoint %>% filter(cjt_safety2 == "Hometown or Syria safe" & cjt_conscription == "Military conscription ended") %>% pull(outcome) %>% table %>% prop.table


#Your hometown is not safe AND military conscription has not ended
basemod_unsafe_conscription1 = lm_robust(outcome ~  cjt_econ_syria_jobs + cjt_econ_lebanon_jobs + cjt_networks, weights = weights, clusters = response_num, data = conjoint %>% filter(cjt_safety2 == "Your hometown is not safe", cjt_conscription == "Military conscription remains")) %>% tidy
basemod_unsafe_conscription2 = lm_robust(outcome ~  cjt_econ_syria_services + cjt_econ_lebanon_services + cjt_networks, weights = weights, clusters = response_num, data = conjoint %>% filter(cjt_safety2 == "Your hometown is not safe", cjt_conscription == "Military conscription remains")) %>% tidy

pevec = c(basemod_unsafe_conscription1[1:2,2], basemod_unsafe_conscription2[2,2], basemod_unsafe_conscription1[5,2], basemod_unsafe_conscription2[5,2], basemod_unsafe_conscription1[8:9,2])
sevec = c(basemod_unsafe_conscription1[1:2,3], basemod_unsafe_conscription2[2,3], basemod_unsafe_conscription1[5,3], basemod_unsafe_conscription2[5,3], basemod_unsafe_conscription1[8:9,3])
pvalvec = c(basemod_unsafe_conscription1[1:2,5], basemod_unsafe_conscription2[2,5], basemod_unsafe_conscription1[5,5], basemod_unsafe_conscription2[5,5], basemod_unsafe_conscription1[8:9,5])

unsafe_reg = lm_robust(outcome ~  cjt_econ_syria_jobs + cjt_econ_lebanon_jobs + cjt_networks, weights = weights, clusters = response_num, data = conjoint %>% filter(cjt_safety2 == "Your hometown is not safe"), cjt_conscription == "Military conscription remains")

p2 = effect_plotter(estimate_vec = pevec, se_vec = sevec, names.variables = c("Jobs in Syria", "Public Services in Syria", "Jobs in Lebanon",
                                                                              "Public services in Lebanon", "Networks"), 
                    names.levels = list(c("Few job opportunities", "Many job opportunities"),
                                        c("Unavailable public services", "Available public services"),
                                        c("Lack good job", "Possess good job"),
                                        c("Public services not available", "Public services available"),
                                        c("Friends/family elsewhere", "Friends/family in Lebanon", "Friends/family in Syria")
                    ),
                    effect.label = "Effect on probability of intending to return", 
                    x.lower = -0.2, x.upper = 0.2, labs = T, title = "Not Safe")
p2


#All Syria is safe AND military conscription has ended
conjoint$cjt_conscription %>% table
conjoint$cjt_safety2 %>% table
basemod_safe_noconscription1 = lm_robust(outcome ~  cjt_econ_syria_jobs + cjt_econ_lebanon_jobs + cjt_networks, weights = weights, clusters = response_num, data = conjoint %>% filter(cjt_safety2 == "Hometown or Syria safe", cjt_conscription == "Military conscription ended")) %>% tidy
basemod_safe_noconscription2 = lm_robust(outcome ~  cjt_econ_syria_services + cjt_econ_lebanon_services + cjt_networks, weights = weights, clusters = response_num, data = conjoint %>% filter(cjt_safety2 == "Hometown or Syria safe", cjt_conscription == "Military conscription ended")) %>% tidy

pevec = c(basemod_safe_noconscription1[1:2,2], basemod_safe_noconscription2[2,2], basemod_safe_noconscription1[5,2], basemod_safe_noconscription2[5,2], basemod_safe_noconscription1[8:9,2])
sevec = c(basemod_safe_noconscription1[1:2,3], basemod_safe_noconscription2[2,3], basemod_safe_noconscription1[5,3], basemod_safe_noconscription2[5,3], basemod_safe_noconscription1[8:9,3])
pvalvec = c(basemod_safe_noconscription1[1:2,5], basemod_safe_noconscription2[2,5], basemod_safe_noconscription1[5,5], basemod_safe_noconscription2[5,5], basemod_safe_noconscription1[8:9,5])

safe_reg = lm_robust(outcome ~  cjt_econ_syria_jobs + cjt_econ_lebanon_jobs + cjt_networks, weights = weights, clusters = response_num, data = conjoint %>% filter(cjt_safety2 == "Hometown or Syria safe"), cjt_conscription == "Military conscription ended")

p1 = effect_plotter(estimate_vec = pevec, se_vec = sevec, names.variables = c("Jobs in Syria", "Public Services in Syria", "Jobs in Lebanon",
                                                                              "Public services in Lebanon", "Networks"), 
                    names.levels = list(c("Few job opportunities", "Many job opportunities"),
                                        c("Unavailable public services", "Available public services"),
                                        c("Lack good job", "Possess good job"),
                                        c("Public services not available", "Public services available"),
                                        c("Friends/family elsewhere", "Friends/family in Lebanon", "Friends/family in Syria")
                    ),
                    effect.label = "Effect on probability of intending to return", 
                    x.lower = -0.2, x.upper = 0.2, labs = F, title = "Safe")

library(ggpubr)
ggarrange(p2, p1, widths = c(1.6, 1)) %>% 
  ggexport(filename = "figures/conjoint_safety2.png", height = 1800, width = 1600, res = 250)


# Export Table ------------------------------------------------------------
texreg(l = list(main_reg, unsafe_reg, safe_reg), digits = 3, custom.coef.names = c("Intercept", "Hometown Safe", "All Syria Safe", "Many Jobs in Syria", "Public Services Unavailable in Syria",
                                                          "Public Services Available in Syria", "Conscription Ended", "Good job in Lebanon", "Public Services Unavailable in Lebanon",
                                                          "Public Services Available in Lebanon", "Friends in Lebanon", "Friends in Syria"), omit.coef = "not available|hard to", 
       override.coef = list(c(basemod1[1:4,2], 0, basemod2[4,2], basemod1[7:8,2], 0, basemod2[8,2], basemod1[11:12,2]),
                            c(basemod_unsafe_conscription1[1:2,2], 0, basemod_unsafe_conscription2[2,2], basemod_unsafe_conscription1[5,2], 0, basemod_unsafe_conscription2[5,2], basemod_unsafe_conscription1[8:9,2]),
                            c(basemod_safe_noconscription1[1:2,2], 0, basemod_safe_noconscription2[2,2], basemod_safe_noconscription1[5,2], 0, basemod_safe_noconscription2[5,2], basemod_safe_noconscription1[8:9,2])), 
       override.se = list(c(basemod1[1:4,3], 0, basemod2[4,3], basemod1[7:8,3], 0, basemod2[8,3], basemod1[11:12,3]),
                          c(basemod_unsafe_conscription1[1:2,3], 0, basemod_unsafe_conscription2[2,3], basemod_unsafe_conscription1[5,3], 0, basemod_unsafe_conscription2[5,3], basemod_unsafe_conscription1[8:9,3]),
                          c(basemod_safe_noconscription1[1:2,3], 0, basemod_safe_noconscription2[2,3], basemod_safe_noconscription1[5,3], 0, basemod_safe_noconscription2[5,3], basemod_safe_noconscription1[8:9,3])), 
       override.pvalues = list(c(basemod1[1:4,5], 0, basemod2[4,5], basemod1[7:8,5], 0, basemod2[8,5], basemod1[11:12,5]),
                               c(basemod_unsafe_conscription1[1:2,5], 0, basemod_unsafe_conscription2[2,5], basemod_unsafe_conscription1[5,5], 0, basemod_unsafe_conscription2[5,5], basemod_unsafe_conscription1[8:9,5]),
                               c(basemod_safe_noconscription1[1:2,5], 0, basemod_safe_noconscription2[2,5], basemod_safe_noconscription1[5,5], 0, basemod_safe_noconscription2[5,5], basemod_safe_noconscription1[8:9,5])), 
       include.ci = F,
       caption = "Model 1: Main conjoint results. Model 2: Conjoint results for vignettes that included hometown is not safe AND military conscription remains. 
       Model 3: Conjoint results for vignettes that included (hometown is safe OR all Syria is safe) AND military conscription ended.",
       label = "tab:conjoint_main",
       file = "tables/conjoint_main.tex",
       single.row = T,
       booktabs = T,
       dcolumn = T,
       caption.above = TRUE,
       use.packages = FALSE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       float.pos = "h!"
)